\chapter{Воксельная графика}

В данной главе рассмотрены различные подходы к компьютерному представлению трехмерных сцен. Так дается обзор технологии CUDA и ее приложений к ускорению операций воксельной графики. Так же рассмотрена полип-ориентированная модель роста коралла, ставшая полем для экспериментов и источником сложной геометрии при выполнении данной работы.

\section{Представление трехмерных сцен}
\subsection{Полигональное}

Полигональное представление трехмерных объектов является на данный момент наиболее распространенным подходом в интерактивной трехмерной графике. При этом подходе, поверхность моделируемого объекта аппроксимируется многоугольными примитивами. Такой подход обладает следующими достоинствами:

\begin{itemize}
\item давно известный и широко применяемый метод, множество публикаций, библиотек и приложений для работы с полигональным представлением;
\item современные графические ускорители способны визуализировать сцены из миллионов треугольников со скоростью, приемлемой для интерактивных приложений;
\item простота геометрических трансформаций и деформации объектов.
\end{itemize}

Однако полигональный подход помимо очевидных достоинств обладает и некоторыми недостатками. В частности, внесение в полигональную модель изменений, изменяющих топологию объекта является достаточно нетривиальной задачей. Так же достаточно нетривиальны в реализации логические операции над полигональными моделями. 

\subsection{Массив вокселей}

Воксельная сцена представляет собой регулярную трехмерную решетку кубических вокселей. Таким образом, воксели можно считать трехмерным расширением пикселей, составляющих двухмерные растровые изображения. Каждый воксель хранит специфичную для приложения характеристику занимаемого им пространства. Таким характеристиками могут быть: цвет, прозрачность, плотность, температура и т.д. Одним из недостатков такого подхода являются очень высокие требования к памяти. Объем памяти, требуемой для хранения кубического объема с разрешением $n$ вокселей на каждое измерение, пропорционален $O(n^3)$, что является существенным недостатком данного метода представления трехмерных объектов.

\subsection{Разреженные структуры данных: воксельное октарное дерево}

Часто при решении некоторой задачи интерес представляет на весь объем воксельной сцены, а только некоторая его область. Например, это может быть граница моделируемого объекта. В этом случае колличество вокелей, необходимых для представления объекта пропорционально $S/h^2$, где $S$ -- площадь поверхности объекта, а $h$ -- линейный размер вокселя. Следует заметить, что исследуемые объекты часто могут иметь фрактальную сруктуру, площадь которой может быть не определена. В этом случае можно рассматривать площадь некой аппроксимации, отбрасывающей детали объекта размером меньше $h$. Площадь такой поверхности можно оценить как $S \in O((l/ h)^D)$, где $l$ -- линейный размер объекта, $D$ -- фрактальная размерность поверхности. С более формальными рассуждениями свойствах фрактальных объектов можно ознакомиться в книге [].
Таким образом, если фрактальная размерность поверхности объекта близка к двум, можно добиться значительного сокращения объема памяти занимаемой сценой, по сравнению с представлением ввиде неупакованного трехмерного массива вокселей.

Для компактного хранения подмножества вокселей сцены требуется некая разреженная структура данных. Примерами таких структура являются:

\begin{itemize}
\item пространственный хэш \cite{gpu_hash},
\item RLE-упакованный массив \cite{voxel_rle},
\item октарное дерево \cite{gpu_octree}, 
\item и т.д.
\end{itemize}

В данной работе рассматривается хранение и визуализация сцены, представленной разреженным воксельным октарным деревом (sparse voxel octree, SVO) \cite{svo_olick}. 

\begin{figure}[h]
  \begin{center}
    \subfigure{}{\includegraphics[height=8cm]{tree-heir.pdf}} 
    \subfigure{}{\includegraphics[height=8cm]{tree-rays.pdf}}
  \end{center}
  \caption{Разрезенное октарное дерево}
  \label{fig:svo}
\end{figure}

Для редактирования сцены, представленной воксельным октарным деревом, в данной работе применяются логические операции над поддеревьями.

\section{Визуализация воксельных объектов}

Существуют различные подходы к визуализации воксльных сцен. В данной работе используется подход, основанный на обратной трассирове лучей.

Идея метода заключается в трассировке от наблюдателя одного или нескольких лучей для каждого пикселя генерируемого изображения. В процессе трассировки отслеживаются столкновения луча с объектами сцены и, на основе их свойств (таких как цвет, прозрачность, освещенность и т.д.)  расчитывается их вклад в финальный цвет пикселя.

\subsection{Объемная визуализация}

В данной работе реализавано два алгоритма визуализации воксельной сцен, работающие с двумя разными из представлениями. Первый из них, реализация которого описана в главе ++, используется для отладочной визуализации поля концентрации планктона для модели просто коралла, описанной в главе ++. Этот метод базируется на методе объемного рендеринга (volume rendering). Суть метода заключается в визуализации некоторого трехмерного скаларного поля $v$ ввиде "облака". Скаларное поле представлено трехмерным массовом, хранящим значения $v$ в узлах регуларной решетки. Для получения значения поля в промежуточных между узнами точках пространства используется трилинейная интеполяция. Для визуализации различных диапазонов значений скалярного поля различными цветами, а так же для выделения только интересующего диапазона используют функицию переноса $T(v) = (q_r, q_g, q_g, \kappa)$, отображающую скалярное поле в интенсивность излучения трех цветовых компонент и коэффициент поглащения.

В процессе рендеринга изображения для определения цвета каждого писеля решается уравнение переноса света вдоль луча, параметризованного переменной $t$: 
$$ I(t) = I(t_0)\tau(t_0, t) + \int_{t_0}^t q(s)\tau(s, t) ds $$
$$ \tau(t_1, t_2) = e^{-\int_{t_1}^{t_2} \kappa(s)ds} $$
где $I(t)=(r, g, b)$ -- интенсивность светового потока в направление наблюдателя в точке луча. Интенсивность расчитывается для трех цветовых компонент $r$, $b$, $g$ -- красного, зеленого и синего. $\tau(t_1, t_2)$ обозначает прозрачность участка луча между точками $t_1$ и $t_2$. $\kappa(t)$ обозначает коффициент поглащения света, а $q(t)$ -- коэффициент излучения.

Для вычисления интенсивности путь луча через сцену разбивается на $n$  сегментов, раздеяемых точками ${t_0, t_1, ..., t_n}$. Количество шагов подбирается так, чтобы были как можно менее заметный артефакты визуализации, возникающие из-за недостаточно точного расчета светового потока, но производительность оставалась приемлимой.

Введем обозначения $$ \theta_k = \tau(t_{k-1}, t_k), \;\;\; b_k = \int_{t_{k-1}}^{t_k}\! q(s)\tau(s, t_k) ds, $$ тогда $$ I(t_n) = I(t_{n-1})\theta_n + b_n = (I(t_{n-2})\theta_{n-1} + b_{n-1})\theta_n + b_n = ...$$ 
$$ = \sum_{k=0}^n b_k \prod_{j=k+1}^n \theta_j, \;\;\; b_0 = I(t_0).$$

Интенсивность луча света, проходящего от внешней границы сцены к наблюдателю можно расчитать при помощи следующего алгоритма:

\begin{algorithmic}
\STATE $intensity = b_0$
\FOR{$k = 1; k \le n; k = k + 1$} 
\STATE $intensity = \theta_k intensity + b_k$
\ENDFOR
\end{algorithmic}

Однако, часто удобнее трассировать луч не к наблюдателю, а от него. Это позволяет, например, отбросить первые члены суммы, вклад которых в результат оказывается меньше некоторого порога, что может существенно улучшить производительность, для сцен, содержащих много непрозрачных участков.

\begin{algorithmic}
\STATE $intensity = b_n$
\STATE $transparency = \theta_n$
\FOR{$k = n-1; k \ge 0; k = k - 1$} 
\STATE $intensity = intensity + b_k transparency$
\STATE $transparency = \theta_k transparency$
\ENDFOR
\end{algorithmic}

Так же обратный порядок обхода удобен при ингерации полигональной и объемной графки.

С деталями реализации данного метода с использованием GPU и результатами работы можно ознакомиться в разделе в второй главе.

\subsection{Освещение воксельной поверхности}

В данной работе рассматривается и другой метод визуализации воксельных объектов при помощи трассировки лучей. Этот метод работает со сценой, в которой все пространство делится на абсолютно прозрачные области между объектами и абсольтно непрозрачные границы объектов. В этом случае для определения цвета пикселя требуется найти только первое пересечение луча, ведущего от наблюдателя, с поверхностью объекта. 

В реализацию данного метода было рещено добавить возможность освещения поверхности точечными источниками света. Использована модель затенения Фонга[], которой, для расчета диффузной компоненты освещения и зеркального блика требуется информация о векторе нормали к поверхности в освещаемой точке. Поэтому каждый воксель поверхности, в использованном представлении сцены, помимо информации о цвете, хранит также нормаль к поверхности. Интенсивность освещенности в данной модели расчитывается как $$ I_p = k_a i_a + \sum_\mathrm{m \; \in \; lights} (k_d (L_m \cdot N) i_d + k_s (R_m \cdot V)^{\alpha}i_s $$
где  $k_a, k_d, k_s$ -- коэффициенты интенсивности амбиентной, диффузной и спектральной компонент освещения для объекта, $k_a, k_d, k_s$, $N$ -- нормаль к поверхности, $L$ -- направление на источник света, $V$ -- направление на наблюдателя, $R$ -- направление отраженного луча, а $\alpha$ -- коэффициент, определяющий резкость бликов.

Так как все пространство мезду границами и внутри объектов игнорируется, для компартного хранения пограничных вокселей используется разреженное октарное воксельное дерево (SVO) ++.

\section{Полип-ориентированная модель роста коралла}
В качестве источника сложной геометрии, а так же для экспериментов с реализацией алгоритмов над воксельными массивами и вокселизацией полигональных сеток, в рамках данной работы была реализована \emph{полип-ориентированниая модель роста коралла}\cite{merks_coral}. 

Пример коралла, сгенерированного при помощи данной модели можно найти на рис.~\ref{fig:coral}. В данной модели поверхность скелета коралла представляется сеткой треугольнгиков. Однако, на каждой итерации алгоритма выполяенся ее вокселизация на регулярную решетку.

\begin{figure}[ht]
\center
\includegraphics[width=50mm]{coral.png}
\caption{Коралл, построенный полип-ориентированной моделью}
\label{fig:coral}
\end{figure}

\section{Особенности вычислительной архитектуры CUDA}

Архитектура CUDA представляет собой набор программных и аппаратных модулей, позволяющих получить значительный прирост производительности при решении многих задач благодаря переносу больших объемов вычислений на графический ускоритель. Программная часть комплекса состоит из набора библиотек, драйвера, а так же SDK, содержащего специализированны компилятор. Аппаратная часть является платой расширения, содержащей 16 мультипроцессоров. Каждый мультипроцессор имеет SIMD архитектуру и способен параллельно исполнять 8 одинаковых инструкций над различными наборами данных, однако с точки зрения разработчика все выглядит так, как будто каждый мультипроцессор параллельно выполняет 32 потока. 

Каждый мультипроцессор имеет небольшое количество быстродействующей памяти на борту, а плате содержится основная оперативная память. Подробнее особенности устройства мультипроцессоров и их взаимодействия с памятью рассмотрены разделе~\ref{c:mem_model}

При программировании при помощи CUDA, GPU рассматривается, как вычислительное устройство, способное параллельно выполнять очень большое число потоков. Ускоритель используется в качестве сопроцессора, на который CPU может передать трудоемкие вычисления, хорошо распараллеливаемые по данным (data-parallel).

Другими словами, функция, которая многократно применяется к независимым порциям данных, может быть выполнена на ускорителе во множестве параллельных потоков. Для этого ее необходимо скомпилирована в набор команд, поддерживаемых вычислительным акселератором, и полученный код, называемый ядром (kernel), должен быть передан на ускоритель.

Центральный процессор и ускоритель работают с различными модулями памяти, называемыми основной памятью (host memory) и памятью устройства (device memory). Для передачи данных между этими видами памяти используется специальное API.

\subsection{Иерархия потоков}
Потоки, выполняющиеся на ускорителе организованы в решетку (grid) блоков (blocks), показанную на рис.~\ref{fig:cuda_grid}.

\begin{figure}[ht]
\center
\includegraphics[width=100mm]{cuda_grid}
\caption{Иерархия потоков}
\label{fig:cuda_grid}
\end{figure}

Блоком потоков называется группа, способная обмениваться данными при помощи быстрой разделяемой памяти (shared memory), а также синхронизировать свое выполнение для координации работы с памятью. Синхронизация производится благодаря указанию в коде ядра специальных точек синхронизации, где выполнение потоков приостанавливается до тех пор, пока все потоки блока не достигнут этой точки. Все потоки одного блока физически выполняются на одном мультипроцессоре, однако каждый мультипроцессоров может выполнять несколько потоков.

Каждый поток имеет собственный идентификатор (thread ID), являющийся номером потока в блоке. Чтобы упростить обработку данных, имеющих многомерную структуру (изображения или объемные массивы данных) имеется возможность задания блока, как двух или трехмерного массива указанного размера. В общем случае потоки обозначаются трехмерными векторами $(x, y, z)$. Идентификатор потока вычисляется по формуле $threadId = x + yD_x + zD_xD_y$, где $D_x$, $D_y$ и $D_z$ обозначают соответствующие размерности блока. Размеры блока задаются в момент запуска ядра и являются одинаковыми для всех блоков решетки.
 
Существует ограничение на количество потоков в одном блоке. Для устройства, использованного при написании данной работы оно равняется 512. Однако блоки одинакового размера, выполняющие одинаковое ядро могут быть объединены в решетку. Таким образом, при вызове ядра можно создать значительно больше потоков. Количество блоков, которое может параллельно обрабатываться акселератором зависит от характеристик конкретного устройства. Если количество блоков превышает параллельные возможности ускорителя, он последовательно начинает обрабатывать новые блоки по мере освобождения вычислительных ресурсов. Подробнее механизм распределения заданий исследуется в разделе \ref{c:cuda_hw}

Так как потоки разных блоков могут выполняться на разных мультипроцессорах, обмен данными между блоками затруднен. Единственным стандартным методом синхронизации потоков из различных блоков, который поддерживался ускорителем, использованным при написании данной работы, было ожидание завершения выполнения всего ядра. Стоит заметить, что в более новых ускорителях реализована поддержка атомарных функций, но здесь они не рассматриваются.

Каждый блок имеет свой идентификатор (block ID). Блоки могут быть организованы в двухмерную решетку указанного размера. В этом случае, блоки обозначаются двумерными векторами $(x, y)$. Идентификатор блока в этом случае вычисляется по формуле $blockId = x + yD_x$, где $D_x$ --- количество столбцов в блоке.

\subsection{Модель памяти}
\label{c:mem_model}
Модель памяти ускорителя изображена на рис.~\ref{fig:cuda_memory}.

\begin{figure}[ht]
\center
\includegraphics[width=100mm]{cuda_memory}
\caption{Модель памяти}
\label{fig:cuda_memory}
\end{figure}

Регистры и локальная память доступны в пределах одного потока, разделяемая память позволяет обмениваться данными потокам, находящимся в одном блоке, а глобальная, константная и память текстур являются общими для всех потоков. При этом константная память и память текстур доступны только для чтения, а остальные виды --- для чтения и записи.

Глобальная, константная и текстурная памяти оптимизированы для различных режимов использования. Дополнительно, текстурная память поддерживает различные режимы адресации и методы фильтрации данных, применяемые в трехмерной графике для работы с текстурами.

При написании вычислительных ядер следует иметь ввиду, опреции с глобальной памятью требуют достаточно много времени, так как доступ к ней не кэшируется. Одним из способов борьбы с этими задержками является способ управление потоками, описанный в разделе \ref{c:cuda_hw} Другим способом является использование быстрой распределяемой памяти в качестве кэша. Другие методы оптимизации описаны в \cite{cuda_struct_alg} и \cite{cuda_optim}.

\subsection{Аппаратная архитектура ускорителя}
\label{c:cuda_hw}

Для написания эффективных програм для CUDA необходимо знать особенности аппаратного устройства ускорителя. Его структура изображена на рис.~\ref{fig:cuda_hw}.

\begin{figure}[ht]
\center
\includegraphics[width=100mm]{cuda_hw}
\caption{Аппаратная структура ускорителя}
\label{fig:cuda_hw}
\end{figure}

Акселератор реализован как набор мультипроцессоров, имеющих SIMD архитектуру. Каждый такт, все процессоры мультипроцессора, выполняют одну и туже инструкцию, но над разными наборами данных. Каждый мультипроцессор содержит четыре вида памяти:

\begin{itemize}
\item Набор 32-х битных регистров;

\item распределенная память;

\item константный кэш, используемый для ускорения доступа к константной памяти, реализованной ввиде региона памяти устройства, доступного только для чтения;

\item текстурный кэш, используемый для ускорения запросов к текстурам, которые хранятся в участке памяти устройства, доступном только для чтения.
\end{itemize}

Глобальное и локальное пространства памяти реализованы как доступные для чтения и записи регионы памяти устройства.

\subsection{Управление потоками}

Группы блоклов (batch) потоков назначаются на исполнение планировщиком заданий. Каждый мультипроцессор последовательно обрабатывает одну группу за другой. Каждый блок обрабатывается только одним мультипроцессором. Благодаря этому разделяемая память размещается в быстродействующей памяти мультипроцессора и работает значительно быстрее глобальной.

Максимальное количество блоков в группе, которую может обработать мультипроцессор, зависит от того, как много регистров использует каждый поток и сколько разделяемой памяти требуется каждому блоку, так как эти ресурсы распределяются между всеми потоками всех блоков, работающих на мультипроцессоре. Если ресурсов недостаточно даже для запуска единственного блока, ядро не может быть запущено.

Блоки, принадлежащие группе, обрабатываемой на мультипроцессоре, называют активными блоками. Потоки активных блоков называются активными потоками. Каждый активный блок разбивается на группы потоков, называемые warp'ами. Каждый warp содержит одно и тоже количество потоков, называемое размером warp'а. Потоки одного warp'а работают в SIMD режиме. Все инструкции этих потоков выполняются синхронно. Если в ядре есть ветвления, и в warp'е есть потоки, идущие в разные ветки, выполнение этих потоков сериализуется. По этой причине расходящиеся ветвления могут снизить производительность.

Warp'ы активных потоков поочередно выполняются мультипроцессором. Ими управляет планировщик потоков, который периодически переключает потоки, минимизируя использование вычислительных ресурсов. Например, когда один warp ожидает данные из памяти или заблокирован, ожидая другие потоки своего блока, потоки других блоков, выполняющихся на данном мультипроцессоре, могут производить математические операции. Это позволяет более рационально использовать вычислительные ресурсы мультипроцессора[].

Количество потоков, арифметические операции которых ускоритель может выполнять параллельно, равняется числу мультипроцессоров умноженному на размер warp'а. Максимальное количество активных потоков равняется кличеству warp'ов на мультипроцессор, умноженному на размер warp'а, умноженному на количество мультипроцессоров. 

На практике такое количество активных потоков достижимо только на достаточно простых ядрах. Требования ядра к регистрам, разделяемой памяти, и размер блока могут ограничить число активных потоков. Отношение реального числа активных потоков к максимально достижимому называется occupancy. Для расчета этой величины по параметрам ядра в CUDA SDK присутствует специальная утилита.

Следует заметить, что для достижения максимальной производительности, при разработке ядер следует знать о большом количестве тонкостей работы ускорителя.

\subsection{Характеристики исследуемой платформы}
При выполнении данной работы тестирование реализации производилось на двух разных ускорителях:
\begin{itemize}
\item NVIDIA GeForce 8800 GTS
\item NVIDIA GeForce GTX 275
\end{itemize}

 Карта NVIDIA GeForce 8800 GTS обладает следующими характеристиками:
 \begin{itemize}
\item 640 мб основной памяти,
\item 14 мультипроцессоров,
\item 32 потока на warp,
\item максимум 24 warp'а на мультипроцессор,
\item максимум 8 блоков на мультипроцессор,
\item 8192 32-х битных регистров на мультипроцессор,
\item 16384 байта разделяемой памяти на мультипроцессор,
\item тактовая частота потоковых процессоров 1.2 Ghz,
\item максимум $14 \cdot 24\cdot 32 = 10752$ активных потоков.
 \end{itemize}

Вторая карта, NVIDIA GeForce GTX 275 обладает следующими характеристиками:
 \begin{itemize}
\item 896 мб основной памяти,
\item 30 мультипроцессоров,
\item 32 потока на warp,
\item максимум 32 warp'а на мультипроцессор,
\item максимум 8 блоков на мультипроцессор,
\item 8192 32-х битных регистров на мультипроцессор,
\item 16384 байта разделяемой памяти на мультипроцессор,
\item тактовая частота потоковых процессоров 1.4 Ghz.
\item максимум $30 \cdot 32\cdot 32 = 30720$ активных потоков.
 \end{itemize}

\section{Использование GPU для ускорения операций с воксельными объектами}

Решение многих задач воксельной графики требует выполнение одной и той же опреции над большим количеством однородных элементов данных. Такие вычисления обычно хорошо параллелятся и ложатся на вычислительную архитектуру графических ускорителей последних поколений.

В данной работе GPU ускорение при помощи технологий CUDA и OpenGL/Cg используется для задач:
\begin{itemize}
\item вокселизаця полигональной сетки на регулярную решетку в модели роста коралла, ++
\item расчет диффузии планктона и его поглащения полипами, ++
\item совместная визуализация полигонального скелета коралла (растеризация) и облака пространственного распределения планктона (трассировка лучей), ++
\item визуализация вокселизированной поверхности, упакованной в октарное дерево, методом обратной трассировки лучей.
\end{itemize}

\subsection{Вокселизация с тспользованием растеризатора GPU}
В рамках реализации модели роста коралла[] был реализован однопроходный алгоритм вокселизации, описанный в работе[]. Этот алгоритм полагается на использование аппаратной растерицазии полигонов. Данная фугкциональность видеокарт не доступна при использовании CUDA, поэтому алгоритм реализован при помощи OpenGL API с использованием языка шейдеров Cg.

Алгоритм получает на вход треугольную сетку. Выходом алгоритма является массив биттов, соответствующих узлам регулярной трехмерной решетки, значения которых которых равняются единице, если узел лежит внутри объема, оганиченного сеткой, и нулю в противном случае. Для корректной работы алгоритма понятия "внутри" и "снаружи" должны иметь смысл, тоесть поверхность не должна иметь разрывов и самопересечений. Параметры сетки, такие как разрешение решетки и координаты вокселизируемого объема задаются перед запуском алгоритма. Разрешение решетки ограничивается объемом памяти и некоторыми характеристиками видеокарты, о которых подет речь ниже.

Для начало кратко опишем используемые при вокселизации возможности современных GPU. На рис.\ref{fig:gl_pipeline} показамы основные стадии обработки данных, выполняемые GPU во время ренеринга изображения.

\begin{figure}[h]
\center
\includegraphics[width=120mm]{gl_pipeline}
\caption{Основные узла конвейера OpenGL API}
\label{fig:gl_pipeline}
\end{figure}

Графический ускоритель проэцирует поступающие вершины на экран при помощи \emph{вершинного шейдера} (\emph{vertex shader}). Затем, на этих вершинах строятся примитивы, которые затем растеризуются в сетку \emph{фрагментов}. В приложениии вокселизации фрагменты считаются эвивалентными пикселям, хотя, строго говоря, это не всегда так. При использовании техники supersampling[] для сокрытия ступенчатого эффекта (aliasing) для каждого пикселя финального изображения может вычисляться несколько фрагментов.

После вычисоения фрагментным шейдером цвета и глубины каждого фрагмента, выполняется отсечение скрытых фрагментов, при помощи z-буфера и смешивание (blending) нового значения для данного фрагмента с содержимым буфера кадра (framebuffer).

Перечислим особенности, имеющиеся в современных GPU, которые позволяют использовать графичский конвейер для эффективной однопроходной вокселизации полигональных сеток.

\begin{itemize}
\item Представление каждого из 4-х (r, g, b, alpha) цветовых каналов изображения в буфере кадра 32-битными целами числами. Таким образом фрагменый шейдер может записть до $32 \cdot 4 = 128$ бит данных в каждый \emph{render target} буфера кадра.
\item Использование до 8 \emph{render target}'ов одновременно. Это позволяет увеличить количество бит информации, возвращаемых фрагментным шейдером до $128 \cdot 8 = 1024$ бит.
\item Возможность использования битовых операций, в частности \texttt{xor}, при блендинге пикселей в буфер кадра.
\item Отключение использования буфера глубины, благодаря чему возможно поподание в буфер кадра фрагментов всех примитивов, не зависомо от их видимости.
\end{itemize}

Используемый метод позволяет производить вокселизацию в решетку, одно измерение которой (назовем его глубиной) не превышает 1024. Разрешение двух других измерений определяется объемом доступной памяти. В модели коралла использовалась решетка размером $256^3$, которая требовала $256^3 = 2^{24}$ бит $ = 2$ мб. Так же проводились экспериметны в решеткой $1024^3$, которая требовала 128 мб памяти видеоускорителя.

Алгоритм воселизации можно описать следующей последовательность шагов.

\begin{enumerate}
\item Привести OpenGL в нужное сотояние:
\begin{itemize}
\item Установить специально подротовленный framebuffer в качестве выходного.
\item Настроить ортогональную проекцию интересующего участка сцены в этот буфер.
\item Отключить использоваие z-буффера.
\item Установить описанный ниже фражментный шейдер.
\item Включить \texttt{xor} в качестве логической опрерации при записи фрагментов в буфер кадра.
\end{itemize}
\item При необходимости установить значения всех битов буфера кадра в 0.
\item Отрисовать все треугольники сцены один раз в произвольном порядке.
\item Использовать полученую вокселизацию, обращаясь к буферу кадра, как к текстуре.
\end{enumerate}

Основная тонкость алгоритма заклютается в используемой фрагментной программе. С ее исходныйм текстом можно ознакомиться в приложении А. Результатом работы этой программы является столбец из 1024 бит, соответствующий одному столбцу объема. Программа построена таким образом, что в возвращаемом столбце все биты, соответствующие узлам решетки, лежащим выше растеризуемого примитива устанослены в 1, а все лежащие ниже равняются 0. Таким образом, операция \texttt{xor} инвертирует состояние всех битов буфера кадра, лежащих выше растеризуемого примитива. После расеризации всех примитивов получается, что биты, соответствующие узлам, лежащие за границей объекта, инвертируются четное число раз, и раняются 0, а биты, соответствующие внитренним узлам решетки инвертируются нечетное число раз, и становятся равными 1.

Данный алгоритм отличается учень эффективным использованием аппаратных ресурсов видеокатры и обеспечивает высокую производительность, делающую его пригодным для растеризации объектов в решетки с высоким разрещением в реальном времени. Авторы утерждают, что время вокселизации объекта, состоящего из примерно 256 тысяч треугольников в решетку с разрешением $1024^3$ занимает порядка 10 мс. Оценки производительности, полученные на схожем оборудовании при выполнении данной работы сравнимы с авторскими, хотя немного уступают им. Возможно, это связано с недостаточно оптимизированной реализацией алгоритма или особенностями методики проведения измерений.

\section{Трассировка лучей в октарном дереве}

Одной из главных частей данной работы является реализация интерактивной обратной трассировки лучей в октарном воксельном дереве. Идея использования октарных деревьев для ускорения трассировки лучей известна уже давно. В данной работе рассматривается реализация алгоримта трассировки при помощи современных GPU и технологии CUDA.

В этом разделе будут описаны идеи, которые легли в основу реализации трассаровщика лучей на GPU. Сама же реализация, с указанием специфичных для конкретной вычислительной архитектуы деталей будет описана в разделе ++.

Пусть весь объем сцены, упаковынный в октарное дерево лежит в единичном кубе, ребра которого выравнены вдоль осей координат. Куб лежит между точками с координатами $(0, 0, 0)$ и $(1, 1, 1)$. Трассируемый луч задается парой векторов $\mathbf{p}_0 = (p_{0x}, p_{0y}, p_{0z})$ -- начальная точка луча и $\mathbf{d} = (d_{x}, d_{y}, d_{z})$ -- его направление.

Определим функцию $\mathbf{p}(t) = \mathbf{p}_0 + t\cdot\mathbf{d}, \; t \ge 0$, отображающую параметр $t$ в точки на луче. Для простоты при обяснении алгоритма трассировки луче в данном разделе будет использоваться двухмерная длоскость и квадродерево вместо октарного. Обобщене на трехмерный случай достаточно тривиально и будет объяснено немного позднее.

Перед началом трассировки обеспечим выполение условия
что $$d_i > \epsilon,\; i \in \{x, y, z\},$$ где $\epsilon$ -- близкая к нулю положительная величина. Для этого произверем следующие операции:

\begin{algorithmic}
\STATE Вход: $\mathbf{p}_0 = (p_{0x}, p_{0y}, p_{0z})$, $\mathbf{d} = (d_{x}, d_{y}, d_{z})$
\STATE $octant\_mask = 0, dim\_bit = 1$
\FOR{$ i \in \{x, y, z\} $}
	\IF {$|d_i| < \epsilon$} 
	    \STATE $d_i$ = copysign($\epsilon$, $d_i$)
	\ENDIF
	\IF {$d_i < 0$} 
	    \STATE $p_{0i} = 1 - p_{0i}$
	    \STATE $d_i  = -d_i$
	    \STATE $octant\_mask\,  |= dim\_bit$
	\ENDIF
 	\STATE $dim\_bit\, *= 2$
\ENDFOR
\end{algorithmic}

\begin{figure}[h]
\center
\includegraphics[width=120mm]{trace_node}
\caption{Определение порядка прохода луча но дочерним узлам ноды октарного дерева.}
\label{fig:trace_node}
\end{figure}
